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Section I 

1.0 Introduction and Summary of Results 

1.1 Introduction 

This report summarizes the work performed under Task B of NASA Contract 
NAS6-2307. The purpose of the work was to characterize the relationship 
between the S-193 Autotatic Gain Control (AGC) data and the magnitude of 
received signal power, which in essence is the ratio of the "peak of the mean 
waveform" to the "mean of the peak (individual) waveform values." In accordance 
with earlier S-193 program documents we use the abbreviation r-factor to describe 
this ratio. The r-factor will be less than unity, and will be a function of 
off-nadir angle, ocean surface roughness, and receiver signal-to-noise ratio 
(SNR). Of these items, the largest variation of r-factor is expected to be 
due to change of angle off-nadir. The effect of ocean surface roughness should 
be much less important and all work reported here is based on a quasi-flat (but 
diffuse scattering) ocean. The digital computer simulation used in the present 
work includes provision for additive receiver (white) noise, but all work to 
date has been for the zero-noise, infinite SNR case. As described in section 
3.0, we specialized to the noise-free case because of the nature of calibration 
data supplied by the S-193 hardware contractor. 

We considered the possibility of conducting either a hybrid simulation 
or an entirely-analog simulation as a way to estimate r-factor, but the diffi- 
culty in either of these approaches is approximately the same; it is difficult 
to adequately shape the range of input process waveforms to be encountered. 

Largely for this reason, we used an entirely digital simulation with the system 
model to be described in the following section. The main body of this report 
discusses the IF and video impulse response functions used, details of the 
input (expectation value) waveforms, and the results to date. An appendix 
provides a more specific discussion of the digital computer programs used. 

1.2 Summary of Results 

During the time in which this r-factor analysis was undertaken, it appeared 
that the Skylab altimeter contractor would provide all calibration data using 
deterministic (non-fluctuating) calibration waveforms and that r-factor would 
explicitly enter into the calibration and data reduction process. In the 
interim we learned that the contractor was Instead attempting to simulate ocean 
scattered signals through: (1) use of noisy triangular waveshapes for the 
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on-tiadlr case, and (2) noisy rectangular approximations to the off-nadir 
waveforms. As a consequence our effort was redirected to yield a comparison 
of these waveforms with theoretically derived and digitally simulated wave- 
shapes. Results to be presented indicate that a correction factor should be 
used with the r-factcr calibration data ob. lined from these approximate 
waveforms. 

There are other complications: Based cn presently available data, the 

S-193 receiver was operating between -7 and 0* C during the SL-2 mission. All 
published calibration data was taken at other temperatures and the available 
data indicates system characteristics change drastically between test tempera- 
tures. A comparison of calibration data taken using deterministic and random 
waveforms (which exists only for the on-nadir case) shows r-factor to be 
about 0.7 at 5* C, and to be near unity at the other two extremes of tempera- 
tures. 

We therefore conclude, based on cur examination of the recent calibration 
data,* that the entire receiver including the age loop, the age attenuator, 
the boxcar detector, and so forth, appears to be grossly nonlinear outside the 
approximate ambient temperature range and for both low and high SNR values 
(although this latter nonlinearity is to be expected to some degree) . We had 
originally expected that the variation of the r-factor with SNR would be an 
important effect; now these other phenomena are seen to provide more serious 
changes. For these reasons, the possibility of calibration data interpolation 
does not appear feasible. Should calibration data not be available near the 
S-193 receiver operating temperature, the only recourse is to record data on 
the back-up hardware, as an assessment of the accuracy of interpolated data. 

The second difficulty is in the shape of the waveforms extracted from 
presently available quick-look data. If further experimental data studies 
(from SL-2, 3, and A) show the recorded waveforms to seriously depart from the 
theoretically computed waveshapes discussed herein, then r-factor should be 
re-computed for these waveshapes using the digital programs documented in the 
Appendix. 

In summary, we recommend that the S-193 calibration data be compensated 
to reflect r-factor perturbations Sue to waveshapes used in the calibration 


*"S-193 Microwave Radiometer /Scat terometer Altimeter," Calibration Data Report, 
Vol. IB, Revision D, Contract NAS9-11195, GE, 22 March 1973. 
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process. We specifically recommend that the corrections shown in Fig. 6 of 
this report be incorporated into S-193 data reduction procedures, subject to 
(1) further examination of measured and computed waiveshapes and r-f actor 
recomputation as required, and (2) acquisition and study of calibration data 
near the actual operating temperature to better understand the Implied 
nonlinear receiver behavior with temperature. 
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Section II 

2.0 Digital Simulation of 'r-f actor Relationship 

2.1 System Ifodel 

Figure 1 (which is intended to be largely self-explanatory) shows the 
overall system model used to estimate r-factors through a digital simulation 
process. The input signal S(t) is the "input waveform" as discussed in a 
subsequent section; S(t) is time varying over times relatively long compared 
to the IF and video response times. For the digital model, the convolutions 
become finite sums of products with the number of terms in the sums chosen 
on the basis of the quantized time increment and on the rate of decay of the 
impulse response function h(t); the choice of number of terms for the present 
case is discussed in the second section. 

The model shows provision for input of both a "signal" S(t) and a noise 
N, and our computer program includes the noise possibility, but all our work 
to date has been done for the case of N « 0 (SNR »*“>). 

By keeping the outputs GRN 3 and GRN^ both zero, and by replacing GRNj^ 
and GRN 2 of Fig. 1 with unit transfer devices (ones for which the output equals 
the input), the model will produce a "deterministic output." This is important 
for veriflng that the system output expectation value has the desired shape. 
There are two ways of obtaining this output expectation value of course; 
either by use of the deterministic model or by running the random model for 
infinitely long times. CleaT.'ly the first is preferable for finite computation 
budgets. The necessary modification to produce the deterministic model merely 
Involves replacement of the Gaussian random number subroutine GAUSS (see 
discussion in the Appendix) by a "fake" GAUSS which returns the input as 
output when called by the rest of the program. 

2.2 Impulse Response Functions h^(t) and hy(t): 

The exact IF and video impulse response functions (h][ and hy) depend upon 
details of the S-19') circuitry which are not readily available; however, the 
central limit theorem assures that if one has a number of individual time 
functions, fj^(t), f 2 (t), ... , fQ(t) which are in general non-negative, 

fi(t) > 0, 

then under fairly general conditions, their convolution will approach a normal 
curve. 






Figure 1. System Model Used for r>factor Simulation Results 
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f(t) - fi(t)*f2<t)* ... *fn(t) 


^ A(0) ^-(t-n)2/2o2 

a-s/SiT 


where 


00 



f(t)dt 


n = + h2 + • • • hn 

CJ “ 0, + 0/| + • . • 0- 

12 n 

2 

and where n^, are the individual mean and variance of f^Ct). Thus for a 
system comprising several stages fj^(t) in cascade, a normal or Gaussian 
impulse response function becomes a progressively better overall system 
description as the nvimber of stages becomes greater. 

This is discussed briefly by A Papoulis,* with an example showing that 
for a cascade of four Identical one-pole systems each having response 


hi(t) = e““*'u(t) 


Hi(u) 


1 

a+jto * 


i - 1, ••• 4. 


the overall system response h(t) 
by 

h(t) = 


‘ hi(t)* -.. *h 4 ^(t) is given approximately 
-(at-4)2/8 

X e 

2V27ct3 


with a system function H(uj) = H]^(uj)x ... xH^Cm) given approximately by 


H(m) « 


1 -j4u)/a -2o)^/a^ 

~4 ® ® 

a 


— ^ 4(o/q 

The e *' is of course the phase shift associated with the h(t) Gaussian 
being centered at tg ■ 4/a. 


*"Systems and Transfoirms with Applications in Optics", New York; McGraw-Hill, 
1968, pp. 78-81. 
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For our problem, we are not too concerned with where the time origin of 
the problem lies, and we will assume an impulse response function of the form 
h(t) ■ with N an amplitude normalization constant to be 

determined later. The time shift to will be chosen merely for convenience, 
and o will be determined from the quoted width of the IF or video filter. 

As described in the Appendix discussing the computer program, the fimctlonal 
form used for the impulse response function is the shifted Gaussian form, 

h(t) = , 

and we need to express this x in terns of a and to express a in terms of the 
3dB width of the video response. The time shifts or phase shifts are of no 
particular importance to the r-factor simulation, and so setting tg •• 0 
temporarily, 

-(at)2/8 -t^/x^ 

e » e 

or T = V*8/a. 

It is easy to show (based on tabulated values for the normal distribution 
function) that this leads to 

4ir fc 

“ “ OaT » 

0.831VT 1 

or X » ; -r— 

47T fg 

where f^, is the system bandwidth to be simulated. For the video bandwidth 
of 5 MHz, this yields 

xy ® 37.36 nanoseconds, 

where the subscript V denotes video (as the subscript I will denote IF) . The 
IF bandwidth is taken as twice the video, so that 

Xj » 18.68 nanoseconds. 

The video response will dominate the system, and as we are going to 
approximate the continuous video impulse response function hy(t) by a finite 
set of weights hy^ spaced a distance 6t apart, 6t for the video weights (and 
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the same $t Is to be used for the IF hj as well as the Input waveform) must 
be chosen r -lall enc.:>{ i to provide several time samples within the Interval 
Ty. For the work reported, we used 6t ■■ 10 nanoseconds. 

The number of weights hj and hy can be chosen by discarding all weights 
of amplitude less than 2% of the greatest weight; the time shifts toy and to^ 
are chosen simply for convenience. Under these various criteria, the followin' 
quantities were used (In addition to the <St, ti, and xy already given): 


to = 80 nanoseconds, ny = 17 


tg = AO nanoseconds, nj = 8. 


Figure 2 shows the video and IF response functions which result from this 
choice of Input parameters. 


2.3 Input Waveforms 

For "he Input waveforms employed, we refer to a February, 1973, Research 
Triangle Institute (RTI) report by Miller, Brown, and Hayne, "Engineering 
Studies Related to Geodetic and Oceanographic Remote Sensing Using Short 
Pulse Techniques," Contract No. NAS6-2135; In particular, see Figures 
on pages 2-46 to 2-50 together with the discussion of these figures. We have 
required that the Input waveforms In this present work reproduce those wave- 
forms of the RTI February, 1973, report (apart from possible amplitude and 
time-origin differences) when the input waveforms are run in the "deterministic" 
program. 

Waveforms shown in this report (hereafter referred to simply as "RTI 
waveforms") include composite IF and video bandwidth effects; therefore, 
direct use of these waveforms in the present simulation program wouxd have 
the effect of including the IF and video characteristics twice. This is 
significant only for the 0® and 0.5® off-nadir waveforms, and those have been 
recalculated with IF and video bandwidth effects removed. For all other 
angles, the waveshape effect of the IF and video is a simple time-shift with 
imperceptible changes in shape. Figure 3 shows the different Input waveforms 
of the present st idy; all have been scaled to the same peak value and plotted 
versus 3 logarithmic time scale. This figure serves primarily to emphasize the 
different running time requirements of different off-nadir angles. 

We chose to Include the trailing edge of each waveform out as far as 1/10 
maximum value. For ^ » 0® this Is achieved by 700 nanoseconds, whereas the 
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Figure 2, Impulse Response Weighting Functions vs. Time 
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C * 8* case requires about 20 nicroseconds to decay enough; thus since DT 
is fixed by hy considerations (with a value of 10 nanoseconds for our work) > 
the 5*8* case will require 20/. 7 “ 29 times the running time of the 5 ■ 0* 
cases . 

Results were also obtained for a triangular waveform (100 nanosec rise - 
500 nanosec fall) and for rectangular waveforms of width 0.4, 1.0, 4.0, and 
14.0 microseconds; these waveforms were i-sed by GE in measuring results for 
the most recent S-193 calibration data report, ar ' our r-factor results for 
RTI waveforms compared to GE triangular or rectangular waveforms will provide 
a set of corrections to the data corrections indicated by tie GE report. 
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3.0 Slroulation Results 

3.1 Discussion 

The previously discussed r<>factor programs have been run for both the 
RTI waveforms and the 6E approximation waveforms (triangular or rectangular) 
for the off-nadir angles 0, O.S> 1.3» 3.0, and 8.0 degrees with at least 250 
simulated pulses averaged for each case. 

Figure 4 shows a typical single pulse result from the simulation for 1.5* 
off-nadir; the time increment DT - 10 nanoseconds used for these computations 
is explicitly shown in the figure. Taking a smaller DT would increase computer 
running times, and the figure suggests that 10 nanoseconds may be a reasonably 
good choice. Figure 5 shows the average of 50 individual pulses for the 1.5* 
case. Notice that the vertical scale is different in the single-pulse and the 
50 pulse figures. The peak of the mean return occurs at 1110 nanoseconds for 
this 1.5* case. 

The results of our runs are presented by Table I. The table also compares 
the r-factors from the RTI waveforms and from the GE approximation waveforms 
(triangular or rectangular) to derive the correction factor to apply co GE*s 
r-factor correction as discussed in the next section. 

Figure 6 shows the same results which are tabulated in Table I. The 
errors indicated (by ± in the table and by bars in the figure) are estimates 
of (plus or minus) one standard deviation of the mean r-f actor as given. Note 
that the correction factor shown in Figure 6 is not an easily describable 
or monotone process. This is because there is no particular relationship 
between the RTI waveforms and the rectangular approximants used by GE. For 
example, the 3dB width of the computed waveforms may match the width of the 
rectangular waveform at one off-nadir angle and not at another. 

Figure 7 reproduces Figure 2B4, pg. 2-58 of the RTI February, 1973, report 
with the present r-factor results added to the figure; this is to make easier 
the comp .ison with the earlier r-factor work. While we present our current 
results as being of an interim nature, and capable of improvement with additional 
running time, this digital simulation method is flexible and powerful enough 
that the present results should be better than any previous computed or simulated 
restilts. 





Table I. Susasary oZ :;«ue II simulation Results 
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Figures. Results of r-factor Computation vs. Of£-nadir Angle 
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Appendix 

Details of Digital Computer Simulation for r-factor Estimation 
Ge-.eral Discussion 

A program listing Is supplied here and a punched source deck Is being 
separately provided for a Fortran IV computer program which Is a straight- 
forward Implementation of the r-factor estimation model supplied In the 
main body of this report. Apart from a subroutine GAUSS, this program 
should be directly transferable to KASA, Wallops Island. We discuss GAUSS 
In the following two paragraphs, and then return to the general discussion 
of the program logic. 

At the heart of the process Is a subroutine GAUSS (IX,S,AM,V) which 
computes a normally-distributed random variable V, of mean AM and standard 
deviation S, at each subroutine call. The integer variable IX must contain 
an odd Integer with nine or fewer digits on the first call to GAUSS: thereafter, 
each call to GAUSS uses the value of IX returned from the immediately previous 
call. The variable IX is included in the output from the simulaticn program 
so that each additional computer run can resume the random number sequence at 
the point at which the last previous run terminated; without the ability to 
start each new run at the last used value of IX, one would In effect be using 
the same sequence of random numbers over and over rather than moving to later 
in the sequence (this sequence length is at least 2^^ t 12 for our computation). 

Subroutine GAUSS uses another subroutine RANDU(IX,IY,YFL) ; RANDU provides 
a floating-point random number YFL which is uniformly distributed on the 
interval (0,1). Every call to GAUSS actually uses RANDU twelve times, approxi- 
mating a Gaussian-dlstrlbuted random variable by the sum of twelve uniformly 
distributed random variables. Source listings for GAUSS and RANDU are also 
provided. Both GAUSS and RANDU have been taken from the IBM Scientific Subroutine 
Package, Version III, and further details and comments are available in the 
publication, IBM Form H20-0205, "System/360 Scientific Subroutine Package 
(360A-CM-03X) Version III Programmer's Manual," (White Plains, New York, IBM 
Technical Publications Department, 1968). Subroutine RANDU Is machine- 
dependent, and Is specific to the word length of the IBM 360 series computer. 

To transfer the r-factor program to NASA, Wallops Island, subroutine GAUSS will 
have to be replaced by the subroutine RAND(AVG,VAR,N,X) which already exists 
at Wallops Island. Since RAND computes an array X(N) of N normally-distributed 
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Figure Al. Logic Flow Diagram for Subroutine XSIMUL, 
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random variables of mean AVG and variance VAR, whereas our program calls 
only a single number at a time, it would probably be appropriate to rewrite 
RAND, renaming it GAUSS and using the same argument list as the present 
GAUSS. This new GAUSS can then use the present Wallops subroutine RANDOM 
for its uniformly-distributed random numbers. 

Returning to discussion of the overall program, the subroutine XSIMUL 
performs the bulk of the work with a part of XSIMUL being written separately 
as XCNEXT. The digital convolution required is performed by XCNVLV; a 
separate entry point XCFILL assists in filling arrays for convolutions. 

The subroutine XAFILL sets up the IF and video impulse response weights, and 
also calls subroutine EPTFIL which provides the input wavefoms. Subroutine 
XREAO provides input data for each run, and the output for results is handled 
by subroutines XHEAD and XOUTPT. 

A general logic flow diagram for XSIMUL is given by Figure A-1. The 
following paragraphs will discuss additional details of subroutines XAFILL 
and EPTFIL, and then the general form of program data input will be stmarlzed. 

Impulse Response Weights Computed by XAFILL 

Subroutine XAFILL computes a set of N weights h^ from the expression 
hi = K exp [ - (^i-^q)^ . ] , i » 1»2, • • • , N 

where ti = (i-l)6t 
N 

and K = 1/ ^ ^ hj, 
i=l 

Thus, hi are derived from a Gaussian function centered at tg and having a 
width T. The normalization factor K is chosen to provide an output equal 
to the input when the input is held at a constant value for a time at least 
as long as N(St. The various relationships in hi(t) are shown on the following 
page. 

Input quantities clearly are t, t, ot, and N. These variables in 
XAFILL are TAU, TO, DT, and N, respectively. Since XAFILL sets up both 
the IF and the video response functions, these are distinguished by the 
labels I and V, and the following correspondence list results: 
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Figure A2. Sketch Showing Relationship Between t, t , 
and h.(t), i=l, 2, . . . ,N. ° 





Quantity 


XAFILL Variables 


i-1, ••• , Nj 

\ 

i=l, , Ny 

s 

ty 

fit 


HI(I), 1=1, ••• , NHI 
TOI 
TAUI 

HV(I), 1=1, ••• , NHV 

TOy 

TAUV 

DT 


The time increment DT is common to both the IF and video characteristics, 
and is supplied from subroutine EPTFIL called by XAFILL; the remaining 
quantities TOI, TOV, TAUI, TAUV, NHI, and NHV are ready by XAFILL as LISTl 
using the NAMELIST procedure available to some Fortran IV compilers- Much 
of our computation has been done using a Teletype ASR-33 as a time-sharing 
terminal, and input via the NAMELIST avoids the awkwardness of formatted 
data input from a teletype. 

Waveform Input From Subroutine EPTFIL 
The purpose of EPTFIL is to fill an array EPT(NE) (labelled EXP(NE) in 
XSIMUL) with discrete values from the input waveform. (Actually because the 
input signal sequence gets squared by the square— law detection, the square- 
root of the input waveform is loaded into the array.) The time increment DT 
was chosen as 10 nanoseconds from consideration of the video response, and 
the number of points to be put into EPT(I) will depend upon the time extent 
of the input waveform's region of Interest. This varies from about 700 
nanoseconds for the on-nadir result to about 30 microseconds for an off-nadir 
angle of S'*; this latter time corresponds to 3000 points in EPT(I). It is 
clearly unnecessary to read in 3000 values from a previously computed waveform, 
and Instead we enter 20 sample values from the input waveform and Interpolate 
between 20 input data point pairs by use of a spline-fit performed by sub- 
routine ESPLIN. 

The Input y vs. t waveform is entered in the arrays TIN(I) and YIN(I), 

1=1, 20. It is necessary that the TIN(I) values be in ascending order, that 
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TIN(l) be less than or equal to DT, and that TIN(20) be greater than or equal 
to (NE X DT) where NE is the total number of points in EPT. YIN and TIN are 
read into EPTFIL through LISTIA as are NE and IMNP. The Integer IMNP is the 
index at which the "peak of the mean" is to be found for the r-factor compu- 
tation. IMNP will be larger than the value of J for which EPT(J) has its 
maximum because of the time-delay of the IF and video characteristics; for 
the IF and video values of our work this index shift was around 11 corres- 
ponding to a time-shift of about 110 nanoseconds. 

To simulate the triangular and rectangular waveforms used by GE for the 
Skylab age calibration data, a simpler replacement EPTFIL subroutine was 
written to provide these simple waveforms directly rather than through use 
of a spline-fit. No listing has been provided for this trivial alternative 
EPTFIL. 

Summary of Data Input to r-factor Program 

The discussion of XAFILL has already pointed' out the use of the Fortran 
NAMELIST feature. Data input to the program consisted first of LISTl 
providing the parameters to calculate impulse response weights HI(NHI) and HV(NHV), 
then second LISTIA setting up the input wavefotm, and then LIST2 to be described 
here. 

LIST2, as read by subroutine XREAD, inputs the program quantities IXR, RN, 

NN, NP, NMN, and MOS. First, NN and MOS are for program features not used in 
this study; NN should be set to zero at the first use of LIST2 and need not be 
entered thereafter, and MOS need not ever be entered. The integer IXR is for 
the random number generator. If no IXR value is entered, the program will 
continue with the current value, starting initially with lXR-4321. A positive 
IXR in will replace the current IXR, a negative IXR in will cause subroutine 
XSIMUL to return to the main calling program (and there to STOP), and a zero 
entered for IXR will cause a jump back to XAFILL (to permit changing HI, HV, 
or the input waveform if desired). Usually, one would wish in the first LIST2 
to explicitly enter the value of IXR last appearing in the printout of the last 
previous program run, and to leave IXR off all subsequent LIST2 data in a 
given run . • 

The value RN allows additive receiver noise to be simulated. All our 
results to date have been for RN»0; once RN is set to zero in LIST2, it need 
not be entered in subsequent LIST2 calls. The remaining quantities NP and 
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NMN are the number of pulses to be averaged for one mean va\Aform, and the 
number of separate mean return waveforms to be run. 




r 


Mi/V' I.'tilh.'d 4 




5NTFT XCFIIL (XI, JXY. TEMP.NH 


KI = nCD (JXY-1 ,NH) 

I? (NF) 3C0, 300,400 

3C^ TFrF(NH)=XI 

3®T0BM 

4Cf' TEPF(NE)=XI 

PHTOFH 

END 


SDERCOTINE XSIMOMSf ,NET,EXi', S£,‘i £YI,TPYy,HI,NHi, IPYV,HV, NHV) 
OTKENSION ET (NET) ,2XP(NE) ,IPYI (NHI) ,TPYu (NUi) ,HI (KHl) ,IFYV (HHY) , 

1 HV(NHV) 


^ ■™ii I 
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0 


to 


r 
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2t 



. 25- 
25' 


27 




29 




30 . 

31 i 
32. 
33- 




3< 

35 


36 



33 


C 


:0 
4! ; 


C 


42 

43 

44 


e 


46 

47 


€ 


51 



r^KFK =26 

IXF=U 321 ; 

3 CALL XAFILL(l:l,NHI,dV,KHV«£ 4 l^«N£,I(!MP) 

100 KXS=IX 8 

WRTTF ( 3 ,S) IXP 

5 ?CFr>T(*''AT STAHT,IXa = ',I 10 //) 

CA.Il XSEAD(IXE, 3 M«N», UPfNHN^.'lOS^MF^NET) 

IF tlXP) 135 , 1 -^ 7,110 

1''-7 TX 5 = KX 3 
GC TC 3 

IFFOT OF IX ?=0 WILL CONTIMOS COaREN'T IXh ASJ JUMP dACK TO 

--TO XAFIIL:TXR .LT. 0 CAU:>ES ivEfUttN FhOa XSIKUL, ndD iXH .GT. 0 
— BNTFRS THE MFH IXR. (DON'T USB IXa=- 37 ) 

1«5 RFTOPN 
11 0 CC NT I NOB 

DC IC^C IHM= 1 ,NHH 

CALI XHEAD(IXE,RN,NK,NP>MaN,flOS) 

SEISES =0 

XPK=C. 

XFK 2 =‘>. 

DC 2 C'' 1 = 1 , NET 
20 r =T(T)= 0 . 

25 <> rtr: 3 gA 

300 TPYVfI)= 0 . 

IFT =1 

BrK= 0 . 

IF (PN-. 0 GCri) 360 , 360,330 
33 n DC 70 '' 1 = 1 , NHI 

CALL GAaSS(IXR,aN,Q. ,?) 

?FYT (T) =V 

CALL GA 0 SS(IXR, 3 H, 0 .,V) 

7 <)n tFYC(I)=V 
JXY =0 

DO fCC 1 = 1 , NHY 
JXY=jrY +1 

IF iFfr.oocTi) 797 ;t 97 ; 7¥8 

797 SI=0, 

sc=''. 

798 CCSTTNDB 

CALL GAOSS(IXR,BN, 0 . «SI) 

CALL GAUSS (I XB,. 3 M ,0 . ,SQ) 

799 CCKTIKOE 

CALL XCNVLV (XI,SI,JX Y ,TPYI,HI, NHI) 

CALL XCNVLV (XQ,SQ,JXY,TPYC# HI, NHI) 

XC=XI*XI+XQ*XO 

80 '' CALL XCFILL(XO,JXY,rPYV,NHV) 

ir'TKN)'T 5 ' 0 ", 45 TJ 7 B 3 U 

850 S 6 = 0 . 

DO SCO 1 = 1 , NN 




i 
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axif=JXY+i 

90? CALI XCWHXTfST.HET.ISl.KPK. SG,aM,IXB« JXY. IPYI.Tt-Yo. H1»KH1. 
1 TFYV,HV,NH¥, KOS) 

GC TC 450 

?60 r>0 1=1. hh i 

?EYI(T)=0. 

40'' TFyq(I)=0. 

jxY=r 

SI? NEXT CCMH2NT 

IF (K*J) 450,45^,405 

U''S PO 4?5 1-1, KN ; 

TF (MCS-1) 41C,410,415 
410 K5T=IPT 

GC TC 430 • 

415 KET=TIT/ROS 

IF <ROD (IET,ROS) ) 435,430,435 

430 FT(KIT)=0. 


430 FT(KIT)=0. 

435 IFT==I?T+1 

THIS SET OF STATBN3HTS KEEPS INCEXING SkAL AS FOR BN .GT. 0 

450 DC 5^0 1-1, NE 

JX1=JXY+1 
SG=5XP (I) 

500 CALL XCNSXTfST.SET.IST.EPK, SG,3N,IX&. JXY.TPYI,TPYj, HI,NHI, 

1 TFYV,HV,KHV, KOS) 

NPL£?S=NPISES+1 

XFy2=XFK2-F3PK*5PK 

XEK=XFK3-EFK 

IF (KFLSES-NP) 250,550,550 

55*^ XHP-NPISES 

EM!«PK=0. 

EFK=XFK/XNP 

V'?FK = XPK2/XNF-EPK*BP K 

OC 6CC 1=1, NET 
6''0 IT (I) =ET (I)/XNP 

EFyFK=ET(IHNP) 

LA7FP POT TAPE F.SCORU PROVISION UERS OS IN XOUTPI 

C.Ul XOOTPT(ET,NET,2XP, ME,EMHPK , EPK, V EPK ,NPLSES, IXR) 

nO'' CCKTINOB 

GO TC 100 


SCEECOTINE XCNEXT (ST, NET, IET,SFK, SG,RN,IXa, JXY, 

1 TPYT,TPYQ, • .,NHI, TFYV.HV.NHV, BCS) 

OTFEKSICN E: NTT) ,Tpyi{NUI) ,IPYW (NHI) ,H1 (N.U) ,TPYV (NHV) , UV (NHV) 


XCNEXT »STEPS ALONG" TO THE NEXT OUTpOT POINT TO 3E COMPUTED 
IF (SG-.OrOOl) 2,2,3 
2 SI=0. 

SC=0, 
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c ! 

^ • 


^ 1 • 

1 

3 i 

3 

SC XC n 

CALL CAUSS(IXR,SG.O. , SI) 


o:: 

6 ; 

a 

5 

CALL GAUSS(IXP,SG,0. ,SQ) 
IF (FN-.QC0‘'1) 10,10,5 
CALL n.AOSS(rXE,RN,0. ,¥) 


7 . 
^ , 


Sr=SI4V 

CALL GAUSS (IXR, SH,0. ,¥) 
SC=SC4V 


^ ; t 

12 


CALL XCMVLV (XI , SI, JX X ,TP XI ,dl , NHI ) 
CALL XCNVLV (Xw,SQ, Jx X,T?Xg, HI , SHI) 
Xr=XI*XT+XC*XQ 


r'^ 

14 

15 ; 

1U 

IF (KCS-1) 14,14,13 

KET=I5T 

GC TC 15 


^ 16i 

C„; 

18 , 

13 

15 

KEl=T®T/nCS 

IF (FCO(IET,HOS) ) 30,15,30 

C.AIL XCNVLV (XV, XO,JXX,TPXV,riV,NHV) 


.-V '91 

''^zol 
2\ . 


IEX=IFT(-1 

ST (KET)=ET(KET)+XV 
IF (XV-SPK) 25,25,20 


olli 

1 

24j 

2? 

25 

3'' 

EFK=XV 

RSXDFN 

CALL XCFILL(XO,JXX,rPXV,NHV) 


• 25t 

• ( 
2 6 ' 

33 

IET=IET+1 

P.EXOBN 

SSC 


. « •; c 

28: C 

30^ _ 


SCBPCUTINS XHEAD(IXa,BN,NN,NP,NMN,MOS) 


‘ r* 

3’’iJ 

^ ,2} C 

w\ 

K«»I1?S GSMS5AL HSADINvJ fOR CASS TO E£ RUH 
WRITE (3,10) IXR 


* 3* 

'W 35i 

36_ 

1? 

F05HAT (1HC/40HC FOLLOWING ASS INPLT PAKAWETSAS (IXR= 110, 1H) ) 

RSC=FN*RN 

H'>ITS (3,20) P.K,RSQ 


j7. 

3si 

39. _ _■ 

23 

3C 

Fr?nA7(1lH N0IS2 FACTOR RN= FI 2 .6 , 1 bH, AND ITS SQUARE- F12.6 ) 
MFITE (3,30) NN,MOS,NF,NHN 

FCPPATC NN=',I3,*, «CS=',I3,', NP=',I4,», AND NaK=*,l4//) 


40 
41 ; 


BETOFN 

FNC 


.. 43 , C 

44 

•is! 


SOBPOOTIHE X0UTPT(Sr,NET,SPT,N2,BENFK,2PK,VSPK,NPLSES,IXR) 
OTEEKSICN ET (NET) ,3PT (NE) 


4S 

r» 

47 - 

48 

XOnTFT HANDLES PRINTOUT OF RESULTS 
OA = E.»fNPK/FPK 


49 

so 

'SI 

IP 

QD=EFK/SNNPK 
WFITE (3,10) IXR 
FOPHAT(5H rXR= 110) 


■ “i 

S3: 

1 

- 54 ' 



1. 


Q 

1 


3 

4* 

5. 

e; 

7} 



22 ! 
«>3 i 


29 


wriT"^ (3,2'') NPLSES,aA 

2?> POFHATf22»OFCLLQHIHG EESUITS FCE 14. 17H R 'TUHHS tlSLC R« F10«b> 
W51TE (3,30) EPK,E»NPK,0B 

3'' ’'OaHIT(5H EPR=S12.6,7H,EHMPK»S12.6,9H,AND 1/fia F10.6) 

SEry=S037 (ABS(V5PKI ) *SIGMf1.,V£gK) 

asIT? (3,35) SEPK 
35 ?€?!»»'»'(• SIGSA(EPK) = *,B12.5) 

W^YTE f3.a0l ME T 

40 FCF«AT(/‘ THE ',I4,' OUTPUT POINTS ARE:*) 

(3,950) 

95 " FCYHA’^‘(//' »♦*« I,;.., OOIPUT. ««♦♦**)/) • 

VI5TTE (3,S00.) (I,BT(I), 1=1, NET) 

933 PCFSRT((* *♦♦•,5(14, • ;',1PE12. 5, •♦♦♦*))) 

W? ITf (3,65) 

65 FCFH>T(1UC/9H ♦/♦/♦/♦ /I HO) 

5ET0BN 

SBC 


SOESCUTIMF XIFILL(HI,MHI,HV,KHV.EET,NE.IHNP) 
CIFEHSION r.I (NHI) ,HV (NHV) ,BPT (HE) 
flAf?FlIST/LIST1/TOI,rOV,TAUI,TAUV, NI,MV 


I 


i 

I 

1 


24l«C 
25* C 
25* ^ 
2-^ 
23 C 

29l 

30j 

31 ! 

♦ 

33i 

14, 

35, 

36_ 

3?r 

3si 

I 

39- 

40i 

i 

4’ ! 

43 ! 

43; 

44 


XAPIIL SETS UP THE IMPULSE BSSPCSSE FUNCTIONS HI(.) FOR THE IP 

AND HV(.) FCR THE VIDEO, AND SETS Of EXPECTATION VALUE 

F05 BAVEFORK DE CALLING EPTFIL. . . 


SFAC (1,LIST1) 

WHITE (3.IIST1) 

NF3=NI 
NHV ,.V 

CAII EPTPIL(CT,SPT,NE.IHNP) 
SUF=0. 

T = P. 

DC 1CO 1=1, NHI 

Tl= (TCI-T)/TAOI 
X = FXF(- (TI*TI)) 

scg=soa+x 

HI (I)=X 
TOP T=74CT 

DC 2C^ 1=1, NHI 

20'' HI (1)=HI (I)/S0H 

sor=o. 


4S- 

46i 

H I 

47j 

43^. 

4»- 

'so. 


T=0« 

DO 3CC 1=1, NHV 
TI=(TCV-T)/TAOV 
X = FXF(" (TI*TD) 


SUM=S0H+X 
HV (I)=X 
3f<0 T=T4DT 


52i 


S3i 

I 

54j 

'■ «6' 

S’ 


J 


I 




3i 


DC lice 1=1, NHV 

ttOn HV (I1=HV (I)/SOB 

DSLETE THIS AFTER DEBUGGING? 

WSITF (3,450) NHI,r0I,TAUI,DT 

U5'' PCFHAT (*f FOLLOWING ABE THE»,I*»,' VALUES I,H(I) FOE IOI=», 

1 TAOI=« ,F6. 3, »,*/* AND CT = *,F6.3//1 

W!>ITF (3,50*') (I,HI(I), 1=1, HHl) 

5'''' FOr»!AT( (• • , • ,F9.5) ) ) 

WFTTF (3,55'') NHV,TO V,TAUV,DT 

55'' FCFHAT ('CFOLLOBING AHE THE*, 14.* VALUES I,H(i) FOE TOV=», 

1 F6.3,», TAUV=»,F6.3,»,*/’ AND DT=*,Fb.3//) 

W9ITF f3,5CC) (I>aV(I), 1=1, NHV) 

DFLETE ABOVE? 

BETOFN ^ 

END 


SUBRCOTINE XREAD (JXR ,XRN,JN, JP,JHK, JOS»N£,N£T) 

NArEIIST/LlST2/IXR,RN,NN,H?,ailN,HCS 

SUBRCOTINE XBEAD HANDLES INPUT OF PABAHETEBS FOB NEXT SIMULATION 
PUN. .. 

IXF=-37 
JOS = 1 
!iCS=0 

READ (1,LIST2) 

IF tIXF+37) 10,20,10 

TXF=-37 SIGNIFIES MO IXR ENTERED ON &LIST2 ; RETAIN PRESENT 

(SEE TEST PEEFCRMED ON IXB IN ROUTINE XSIMOL...) 

JXE = IXR 
20 XFN=FN 
JN=NK 
JF = NF 
JRK=NNM 

IF (PCS) 40,40,30 
31 JCS=HO“S 
40 NET= (JN+NE)/JOS 
PFTOFN 
END 


SOEPCDTINE EPTFIL(Dr,YOUT,KE,KMNP) 

DTPSNSIOH yOUT(KE) , TOUT (3000) ,YIN (20) ,TlR (20) (5C) ,X (' 

NAFEIIST/LIST1A/IMNP,NE,TIN, YIN 
,.FOP 8 DFGBSES OFF NADIR, MAIN PROG DIMENSIONS UP TO 3000 


56 

57 







C3 r 

I 


o 


31 


V 


• >: C 


4 iC 


SOayOOlINE EPTFIL ossa a SPLINS PaUCTION pit to 20 input WAVEFOEa 
PAI'^S TIK (I) ,YIN (I) ;TH£SE SAtlPLE POINIS CAE BE CHOSEN IC 
C PDECI'AISLY P.EPPESSNT THE ENTIRE INPUT HAVEFOoH. 

C 

DC 15 1 = 1 , KS 
15 YCOT(T)= 0 . 

DT= 10 . 


o 


o 


1 c 

o,. 


TT= 1 C. 

PnPC ( 1 ,LIST 1 A) 
KJ'hP = IMNP 


o!!! 


•S^ 


DC 99 1 = 1,20 
X(I)=TIN(I) 
09 YUl-YINfl) 


, 6 i KE = NE 

0,7! IF (KF- 30 CC) 19 , 19,16 

If NFTTF ( 3 , 17 ) NE 


o 


17 FCBR»T(//* SEROai N£=',I 4 ,* BUT 3000 IS HAXIHUttl wQIT|i') 
STCP 

7 i : 19 DC 2 C 1 = 1 ,N 3 ; 


19 

20 j 


22. TCCT (T) =TI 

©23* 2f- TT=TT4DT 

2^i CALL ESPLIN fX , Y.20 , TOUT , YCUI , NE, « OCO0 1) 


25! 

2 si 

28; 

29 l 

30 ^ 

3 <! 


DC 25 1 = 1 , NS 
XQ = YCUT (I) 

IF (XCI 3 ", 25,_35 


30 YCUT|I)= 0 . 

GC TC 25 

35 YCUT (T 1 =SQ?T (XO) 


25 CONTINUE 


©,2!C TCOl NCW A FREE ARRAY, USB 10 OUTPUT THE INPUT SiyUENCE 

33] DC 1 CD 1 = 1 . NE 

* 34' 

O 3S‘ 


36L 


XC = YCnT(I) 

100 TCOT<I)=XQ*XQ 
WPITE ( 3 , 15 0 ) 


3,i ‘* 5 "' FCFBAT(/» FOLLOWING IS INPUT SEQUENCE, TIMS= 1*(10 NANOSEC)*//) 

Oasi WFITE ( 3 , 200 ) 

2 QD FCPHATr* ****. 5 f*.. I.;, ..OUTPUT, ,.**♦*)/) 


40 - WFITE ( 3 , 25 ^) (I, TOUT (I) , 1 = 1 , NE) 

Onii 250 F 0 rtMAT((' ,5 (14 , • ; • , lESI^:. 5 , » ♦*♦• ) ) ) 


I 

42l 


9 ETDFN 


O44. C 

O 471 


END 


48: 


SUBEOUTINE ES PLI N (X, Y , N, T ,SS , M , EPSLN) 

DTPEKSICN T(M) ,SS(«) 

DIPEW5I0N X( 50 ) ,Y( 50 ) ,H^ 50 ) ,D£LY ( 50 ) ,H 2 ( 50 ) ,B(SG) ,DELSUY ( *' ) , 


^ 49i 

Oso' C 


o; 


so 

51 j ^ 

S2i 


1 52 ( 5 -') ,C( 50 ) ,S 3 ( 50 ) 

THIS IS A SPLINE-PIT ROUTINE AD"PTE£ FROM "SPLINE FUNCTIONS, 


^ 54' 

.ssf 







o 


©::ic 


0„i c 


c ^ 

©..i 


INTFPPCLATICN, AND NUaKfilCAL iiU ADBAIOaE," I. M. E.GaEVILLi, IK 
••lATKEPATICAl KETHOD3 FCB DIGITAL CCttPOTLRS. VOL.II.” A.ivALSiON 6 


H.S.SI1F,ED*S. ,(KEH YOBK; J. BILEY, 1967) , PP. 156-168 

IF (K-50) 2,2,999 


>9 BFITE (3,998) H 

>8 FCFHAT(//» ERRORI HAD‘,14,* INPOT POINTS, E3PLIN ilANULES CNLY 50*/) 
RETOBN 


2 N1*N-1 

3 DC 51 1=1, N1 
H(T)=X(I+1)-X(I 


DEIY (I) = (Y(I+1)-Y(I))/H(I) 
0 DC 52 1=2, N1 

H2(I)=H(I-1)+H(I) 


3 (I) = .5*F (I-1)/H2(I) 

DEISCY (I) = (DELY (I) -DELY (1-1) ) /H2 (I) 
S2 (I)=2.*DEL£QY(I 


>2 C(1)=3.*DELSCY(I) 
S2 (1)=0. 

S2 (N)=0. 


0MFGA=1. 071797 
MI7B=1 

5 ETA=C« 

6 DC 10 1=2, Hi 

7 H= (C (I) -B (I) *S2(I-1) - (.5-fi(I) ) *S2 (1+1) -S2 (I) ) ♦OttBGA 

8 IF {AES(Wl-ETA) 10,10,9 

9 ETA=ABS(H) 

ID S2 (I)=S2 (I)+H 

13 IF (FTA-EP3T.K) 14,14,997 

)7 IF (NIT8-10) 996,996,995 

)6 NITF=NITP+1 

GC TC 5 

15 WFITE (3,994) NITR,BTA 

>4 FCFHAT(/* NUMBER OF IIEEATIONS= • , 13 , • ANO EPSILON= • , E 12 . 5//) 

14 DO 53 1 = 1, Hi 

i3 S3 (I) = (S2 (I+1)-S2(I) )/H (I) 

COEFFICIENTS NOW ESTABLISHED, POLLOHING STEPS FERFOEM DESIRED 
INItRMECIA'fE PCINT INTEfiPOLATION. 

15 DC 61 J = 1,M 

rrr^T 

>4 IF (T(J)-X(I)) 58,17,55 
>5 IF (T(J)-X(N)) 57,59,58 
>6 I F (T ( J ) - X ( I ) ) 6 0 , 1 / , 5 7 
)7 1 = 1+1 
GC TC 56 

)6 WRITE (3,46) J 

(4 FCFHAT(//* *,I4,*TH ARGUMENT OUT CF' RANGE* ,//) 

SS(J)*C. 





\ 


2 GC TC 61 

3 59 I=!i 

« 60 1*1-1 

s 17 HT1*T (J)-X(I) 

6 - HT2=1 

7 PFCD*HT1*HT2 . 

8 SS2*£2(I)+HT1*S3(I) 

9 PllSC5=fS2(IHS2 fH-n 4SS2)/6. 

,0 SS (J)=Y(I)+H11*DELX(I) + PB0D*DELSQS 

i,‘ 61 CCKTINDE 

,2^ BE7UF.N 

I __ 

1 


I 

«6| 





SOPROCTIHS GAUSS (IX, S, AH, V) 



DO 5^' 1*1,12 
CAll BAHDO(IX,IY,Y) 

IX=TY 

5'' A*A+Y • 

V* (A-6.0) ♦S+AH 

R?»rURH 

SOBRCOTINE 5 AMDU (IX, I Y,YFL) 

TT=IX*6S’S39 

TF (lY) 5,6,6 

5 TY=TY42147tt83647+1 

6 YFI*IY 

YFL=YFL*.ti6566l3E-9 

RETDPK 

EPC 



